function code=code_final_submission(CE)

%% Exogenous parameters set in the run file
global eS eB gammaR1 gammaB A3g omega3g delta xihigh xilow cf1 cf2 cf3 wE wB wS ...
     rho  mc  ceq  phiI cI CC phiD cud cud2 wN Rout...
     CR LCR LevR NSFR LR LIQreg scale InCR InLCR InLevR InNSFR InLR InLIQreg InScale...
     Dscale LIQscale Iscale Escale IntD IntLIQ IntI IntE InLIQS

%% Define endgoenous variables
I=CE(1); %loan to entrepreneur
L=CE(2); %liquid asset holdings at t=1
D=CE(3); %deposits
E=CE(4); %equity
xistar=CE(5); %run threshold
rI=CE(6); %loan rate (net). Gross rate in the paper is R_I=1+rI
r3D=CE(7); %late deposit interest rate (net). Gross rate in the paper is R_D=1+r3D
r2D=CE(8); %early deposit interest rate (net). Gross rate in the paper is r_D=1+r2D
psiBS=CE(9); %L multiplier on balance sheet constraint
psiGG=CE(10); %L multiplier on global game constraint
psiLD=CE(11); %L multiplier on loan demand schedule
psiDS=CE(12); %L multiplier on deposit supply schedule
if InLIQS==1
LIQ1R=0; %Savers' liquid asset holdings    
else
    LIQ1R=CE(13);
end
slackSP=CE(14)^2; %Slack participation constraint of banker--for planner's problem
psiSP=CE(15); %L multiplier on participation constraint of banker

% Additonal varaibles for regulatory equilibrium
lambdaCR=CE(16); %L multiplier on capital requirement
lambdaLCR=CE(17); %L multiplier on liquidity coverage ratio requirement
lambdaLevR=CE(18); %L multiplier on leverage ratio requirement
lambdaNSFR=CE(19); %L multiplier on net stable funding ratio requirement
lambdaLR=CE(20); %L multiplier on liquidity to loans ratio requirement
lambdaLIQreg=CE(21); %L multiplier on liquidity to assets ratio requirement
tscale=CE(22); %Pigouvian tax on scale implemented as a tax on deposit interest payment
tD=CE(23); %Pigouvian tax on deposits
tLIQ=CE(24); %Pigouvian tax on liquid asset holdings
tI=CE(25); %Pigouvian tax on loans
tE=CE(26); %Pigouvian tax on equity

%% Useful Auxilary variables
q=(xistar-xilow)/(xihigh-xilow); %Definition of run probability

%% Preference for liquidity parameteters
if phiD==1
gD=cud2+(log(cud+D*(1+r2D)))*ceq;
else
    gD=cud2+((cud+D*(1+r2D))^(1-phiD)/(1-phiD))*ceq;
end
dgD=(1+r2D)*(cud+D*(1+r2D))^(-phiD)*ceq;
d2gD=-phiD*(1+r2D)^2*(cud+D*(1+r2D))^(-phiD-1)*ceq;

%% Utilities

% % Banker's utility function
UB=rho*U1(eB-E,gammaB,cf1)...
  +integral(@(xi)omega3g*((xi*I-delta*D*(1+r2D)+L)/xi*(1+rI)-(1-delta)*D*(1+r3D))-mc,xistar,xihigh,'ArrayValued',true)/(xihigh-xilow);
% % Derivatives of UB wrt:
UBdI=omega3g*(1+rI)*(1-q);
UBdD=-omega3g*(delta*(1+rI)*(1+r2D)*log(xihigh/xistar)/(xihigh-xilow)+(1-delta)*(1+r3D)*(1+tscale)*(1-q));
UBdLIQ=omega3g*(1+rI)*log(xihigh/xistar)/(xihigh-xilow);
UBdE=-rho*MU1(eB-E,gammaB,cf1);
UBdxistar=-(omega3g*((xistar*I-delta*D*(1+r2D)+L)/xistar*(1+rI)-(1-delta)*D*(1+r3D))-mc)/(xihigh-xilow);
UBdRI=omega3g*(I*(1-q)-(delta*D*(1+r2D)-L)*log(xihigh/xistar)/(xihigh-xilow));
UBdRD=-omega3g*(1-delta)*D*(1-q)*(1+tscale);

% % Banker's investment and utility in autarky
Iaut=max(eB-(rho/Rout)^(1/gammaB),0);
UBaut=rho*U1(eB-Iaut,gammaB,cf1)+Iaut*Rout;
    
% UEstar=(xihigh-xistar)/(xihigh-xilow)*(cI*(phiI-1)*I^phiI);
WEdI=(xihigh-xistar)/(xihigh-xilow)*(cI*phiI*(phiI-1)*I^(phiI-1));
WEdxistar=-1/(xihigh-xilow)*(cI*(phiI-1)*I^phiI);

% USstar=U1(eS-D-LIQ1R,gammaR1,cf1)+MU1(eS-D-LIQ1R,gammaR1,cf1)*(D+LIQ1R)+(xihigh-xistar)/(xihigh-xilow)*(gD-dgD*D*(1+r2D));
WSdD=-M2U1(eS-D-LIQ1R,gammaR1,cf1)*(D+LIQ1R)*InLIQS-d2gD*D*(1+r2D)^2*(xihigh-xistar)/(xihigh-xilow);
WSdxistar=-1/(xihigh-xilow)*(gD-dgD*D*(1+r2D));

%% Global game function

%%%Preliminaries: thetastar, hatlambda (hatl), and its derivatives %%%%%%%% 
thetastar=(L+xistar*I)/(D*(1+r2D));

hatl=((xistar*I+L)*(1+rI)-xistar*(D*(1+r3D)+mc/omega3g))/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)));

hatldI=xistar*(1+rI)/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)));

hatldD=-xistar*(1+r3D)/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)))...
      -hatl/D;

hatldLIQ=(1+rI)/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)));

hatldxistar=(I*(1+rI)-(D*(1+r3D)+mc/omega3g))/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)))...
           +hatl*(1+r3D)/((1+r2D)*(1+rI)-xistar*(1+r3D));

hatldRI=((xistar*I+L))/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)))...
       -hatl*(1+r2D)/((1+r2D)*(1+rI)-xistar*(1+r3D));

hatldRD=(-xistar*D)/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)))...
       +hatl*xistar/((1+r2D)*(1+rI)-xistar*(1+r3D));

% hatldr2D=-hatl*D*(1+rI)/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)));

% % % % % GG function and its derivatives % % % % % % % 
                 
GG=omega3g*D*(1+r3D)*(hatl-delta)...
  -D*(1+r2D)*(thetastar-delta)...
  -integral(@(m)(xistar*I+L)/m,thetastar,1,'ArrayValued',true);

GGdI=hatldI*omega3g*D*(1+r3D)...
    -integral(@(m)(xistar/(m)),thetastar,1,'ArrayValued',true);

GGdD=omega3g*(1+r3D)*(hatl-delta)...
    +hatldD*omega3g*D*(1+r3D)...
    -(1+r2D)*(thetastar-delta);  

GGdLIQ=hatldLIQ*omega3g*D*(1+r3D)...
      -integral(@(m)(1/(m)),thetastar,1,'ArrayValued',true);
    
GGdxistar=hatldxistar*omega3g*D*(1+r3D)...
        -integral(@(m)(I/(m)),thetastar,1,'ArrayValued',true);
   
GGdRI=hatldRI*omega3g*D*(1+r3D);
      
GGdRD=omega3g*D*(hatl-delta)...
     +hatldRD*omega3g*D*(1+r3D);
         

%% LD: Derivatives of the Loand Demand

LDdI=omega3g*(A3g-(1+rI))*((delta*D*(1+r2D)-L)/(I^2)*log(xihigh/xistar)/(xihigh-xilow))-(1-q)*cI*phiI*(phiI-1)*I^(phiI-2);

LDdD=-omega3g*(A3g-(1+rI))*(delta*(1+r2D)/I*log(xihigh/xistar)/(xihigh-xilow));

LDdLIQ=omega3g*(A3g-(1+rI))*(1/I*log(xihigh/xistar)/(xihigh-xilow));

LDdxistar=omega3g*(A3g-(1+rI))*(-1/(xihigh-xilow)+(delta*D*(1+r2D)-L)/(xistar*I)/(xihigh-xilow))+1/(xihigh-xilow)*cI*phiI*I^(phiI-1);

LDdRI=-omega3g*((1-q)-(delta*D*(1+r2D)-L)/I*log(xihigh/xistar)/(xihigh-xilow));


%% DS: Derivatives of the Deposit Supply

DSdI=q*(cf2*delta+cf3*(1-delta))*(xistar+xilow)/2/D;           

DSdD=InLIQS*M2U1(eS-D-LIQ1R,gammaR1,cf1)...
    -q*(cf2*delta+cf3*(1-delta))*(L/D+I/D*(xistar+xilow)/2)/D...
    +(1-q)*d2gD;                         

DSdLIQ=q*(cf2*delta+cf3*(1-delta))/D;                         

DSdxistar=1/(xihigh-xilow)*...
            ((delta*cf2+(1-delta)*cf3)*(L/D+I/D*xistar)...
            -delta*cf2*(1+r2D)...
            -(1-delta)*cf3*omega3g*(1+r3D)...
            -dgD);                    

DSdRD=(1-q)*cf3*omega3g*(1-delta);
                                   

%% Equilibrium Conditions

% % Foc wrt I 
code(1)=(wB+psiSP)*UBdI-psiBS+psiGG*GGdI+psiLD*LDdI+psiDS*DSdI*CC...
       +wE*WEdI...
       -lambdaCR*CR+lambdaLCR*xilow-lambdaNSFR*NSFR-lambdaLR*LR-lambdaLIQreg*LIQreg-tI;
       
% % Foc wrt D 
code(2)=(wB+psiSP)*UBdD+psiBS+psiGG*GGdD+psiLD*LDdD+psiDS*DSdD...
       +wS*WSdD...
       -lambdaLCR*LCR*(1+r2D)-lambdaLevR*LevR+lambdaNSFR*wN*(1-delta)-tD;

% % Foc wrt L 
if wE==0.2 && wS==0 
code(3)=L; % % setting L=0 when code yields L<0 for very high wE in SP
else 
    code(3)=(wB+psiSP)*UBdLIQ-psiBS+psiGG*GGdLIQ+psiLD*LDdLIQ+psiDS*DSdLIQ*CC...
       +lambdaLCR+lambdaLR+lambdaLIQreg*(1-LIQreg)-tLIQ;
end

% % Foc wrt E
code(4)=(wB+psiSP)*UBdE+psiBS...
       +lambdaCR+lambdaLevR*(1-LevR)+lambdaNSFR-tE;

% % Foc wrt xistar
code(5)=(wB+psiSP)*UBdxistar+psiGG*GGdxistar+psiLD*LDdxistar+psiDS*DSdxistar*CC...
       +wE*WEdxistar...
       +wS*WSdxistar;
                                 
% % Foc wrt to RI
code(6)=(wB+psiSP)*UBdRI+psiGG*GGdRI+psiLD*LDdRI;                    

% % Foc wrt to RD
code(7)=(wB+psiSP)*UBdRD+psiGG*GGdRD+psiDS*DSdRD;
              
% % Foc wrt r2D
code(8)=r2D;

% % Balance Sheet constraint
code(9)=I+L-D-E;

% % Global Game constraint
code(10)=GG;

% % Loan Demand
code(11)=omega3g*(A3g-(1+rI))*((1-q)-(delta*D*(1+r2D)-L)/I*log(xihigh/xistar)/(xihigh-xilow))-(1-q)*cI*phiI*I^(phiI-1);

% % Deposit Supply
code(12)=-MU1(eS-D-LIQ1R,gammaR1,cf1)...
        +q*(cf2*delta+cf3*(1-delta))*(L/D+I/D*(xistar+xilow)/2)...
        +(1-q)*(cf2*delta*(1+r2D)+cf3*(1-delta)*omega3g*(1+r3D))...
        +(1-q)*dgD;

% % S self-insurance, LIQ_1^S---CE(13) can be negative implying LIQ_1^S=0
 code(13)=-MU1(eS-D-CE(13),gammaR1,cf1)+(cf2*delta+cf3*(1-delta));
             
%% Social planner's contraint: UB>UBaut
code(14)=UB-UBaut-slackSP; 

code(15)=psiSP*slackSP;


%% Regulation
if InCR==1
code(16)=E-I*CR;
elseif InCR==0
code(16)=lambdaCR;
end

if InLCR==1
code(17)=L+xilow*I-LCR*D*(1+r2D);
elseif InLCR==0
code(17)=lambdaLCR;
end

if InLevR==1
code(18)=E-(E+D)*LevR;
elseif InLevR==0
code(18)=lambdaLevR;
end

if InNSFR==1
code(19)=(E+wN*(1-delta)*D)-NSFR*I;
elseif InNSFR==0
code(19)=lambdaNSFR;
end

if InLR==1
code(20)=L-I*LR;
elseif InLR==0
code(20)=lambdaLR;
end

if InLIQreg==1
code(21)=L-(I+L)*LIQreg;
elseif InLIQreg==0
code(21)=lambdaLIQreg;
end

%% Pigouvian taxation
if InScale==1
code(22)=I-E-scale;
elseif InScale==0
code(22)=tscale;
end

if IntD==1
code(23)=D-Dscale;
elseif IntD==2
code(23)=tD+0.022;
elseif IntD==0
code(23)=tD;
end

if IntLIQ==1
code(24)=L-LIQscale;
elseif IntLIQ==0
code(24)=tLIQ;
end

if IntI==1
code(25)=I-Iscale;
elseif IntI==2
code(25)=tI+0.0433;
elseif IntI==0
code(25)=tI;
end

if IntE==1
code(26)=E-Escale;
elseif IntE==0
code(26)=tE;
end


end

